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ABSTRACT 


POWER SPECTRA OF GEOID UNDt LATIONS 


Data from spacecraft altimeters are expected to contribute to an improved de- 
termination of the marine geoid. The geoid is typically described in terms of 
undulations, i.e. t the distance between the geoid and a reference surface, and these 
undulations may be recovered from altitude data and knowledge of the spacecraft 
orbit. To better define altimeter system design requirements for geoid recovery, 
amplitudes of geoid undulations at short wavelengths are examined. Models of 
detailed geoids in selected areas around the earth, developed from a combination of 
satellite derived spherical harmonics and l°-by-l° area mean free-air gravity anom- 
alies, are subjected to a spectral analysis. The resulting undulation power spectra are 
compared to existing estimates for the magnitude of geoid undulations at short wave- 
lengths. The undulation spectra are found to be consistent with Kaula's rule of thumb, 
following an inverse third power relationship with spatial frequency for wavelengths at 
least as small as 300 km. The requirements imposed by this relationship on altimeter 
accuracy, data rate, and horizontal resolution to meet the goal of a detailed geoid 
description are discussed. 
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POWER SPECTRA OF GEO ID UNDULATIONS FOR THE NORTH ATLANTIC 


INTRODUCTION 

Data from spacecraft altimeters are expected to contribute to an improved determina- 
tion of the size and shape of the marine geoid. If the orbit of the spacecraft is known, 
then altitude data provides a direct estimate of the geocentric radius to the geoid 
(References 1,2). Even assuming realistic uncertainties in knowledge of the orbit 
and geopotential parameters, it is theoretically feasible to recover geoid information 
in a bootstrap process (References 3, 4, 5, 6). 

The shape of the geoid is described in terms of geoid undulation, the normal distance 
between the geoid surface and some reference surface (Reference 7). The resolution 
with which this surface can be defined by altimeter measurement is dependent on 
altimeter design and operation. Three factors related to altimeter hardware and data 
processing impose practical limits on the detection of small geoid features. First, 
the geoid undulation must be of sufficient amplitude to be detected by the altimeter; i.e., 
the measurement error of the altimeter must be less than the geoid undulation ampli- 
tude. Second, referring to Figure 1, the effective illuminated surface area (the area 
illuminated when the detected wavefront intersects the surface) must be smaller than 
the wavelength of the undulation to be measured. Finally, the sampling rate of the 
altimeter, l/r g (sec - *), must be such that two measurements can be made within one 
wavelength of the undulation to be measured (Reference 8). For example, the feature 
depicted at the far left of Figure 1 probably would not be resolved by the altimeter due 
to effective illuminated area constraints as well as data rate constraints. 

To evaluate these factors in a simulated operational environment, the best avail- 
able geoid models, i.e., detailed gravimetric geoids developed by Vincent, et al 
(Reference 9 10), were chosen for undulation power spectra analysis. Several authors, 
notably Kr.ula (References 11, 12) and Rapp (References 13, 14) have analyzed the power 
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spectra of gravity anomalies and the gravity potential from observed and modeled 
gravity data. These gravity spectra imply a certain shape for the undulation power 
spectrum on a world-wide basis (Reference 15). However, a smaller quantity of less 
detailed gravity data was used in these analyses than is incorporated in the geoid 
models selected for this analysis. 

A capability for computation of autocovariance of geoid heights in two dimensions which 
eliminates the need for averaging and fitting of curves to one -dimensional profiles, 
was used in the spectral analysis of nine local geoid areas. These spectra are com- 
pared to the spectra of related geophysical parameters, and appropriate conclusions 
are drawn for modeling the spectrum and for altimeter design and operation. 

ANALYSIS 

Detailed Geoid Model 

The geographical areas chosen for evaluation cover about one-fourth of the globe. These 
areas, depicted in Figure 2, are described as the North-East Pacific, Canada, United 
States, South America, North Atlantic, Europe, Eurasis, India, and Australia. For 
these areas geoid heights were available at 1° intervals. These heights were generated 
by combining l°-by-l° area mean free-air gravity anomalies with gravity anomalies 
calculated from the spherical harmonic coefficients of the SAO'69 geopotential model 
and using Stokes' formula following the method used by Vincent, et al (Reference 9). 

The geoid undulation N is expressed by Stokes' formula as 

2 IT TT/2 

N(0,X) = ~^ J J Ag(0',X')S(*) cos 0' d 0' d X' (meters) (1) 
A'=0 0'=-tT/2 

where 0,X are the latitude and longitude, respectively, of the computational 
point in radians 

0\X' are the latitude and longitude, respectively, of the variable integra- 
tion point in radians 
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R 

y 

Ag(0’,X») 

S (*) 

♦ 

This integral is carried out over the entire spherical surface defined by the mean 

radius R. Because of the behavior of Stokes' function, however, the effect of gravity 

anomalies in regions distant from the evaluation point (0,X) is negligible. Accordingly, 

the region of integration is divided in equation (1) into two areas; a local area 

surrounding (0,X) and the remainder of the sphere, denoted A . The anomalous gravity 

z 

field in each region is also divided into two parts; a gravity anomaly Ag corresponding 

s 

to low order and degree spherical harmonies derived from satellite tracking data, and 

a gravity anomaly Ag defined by 

z 

A? (0,X) = Ag (0,\) + Ag (©, X) (2) 

S 2 

where Ag is the mean free-air gravity anomaly at (0,X). 

Thus equation (1) can be rewritten as: 

N(<£, X) = N (0, \) + N o (0, \) + N (0 , \) < 3 ) 

* o 


is mean earth radius (m) 

is mean gravitational acceleration at the earth’s surface in 
m/ sec^ 

2 

is the mean free-air gravity anomaly at (0', X’) in m/sec 



is the great circle central angle between (0, X) and (0',X') in radians 
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SHADED REGIONS ARE GEOID MODEL AREAS 
FOR GEOID HEIGHT SPECTRAL ANALYSIS 


AREAS REPRESENTED BY DETAILED GEOID MODELS 





WtW9W''/ */fA 

















1 y. 

n 

1 // 

Bit 






where 


2tt tt/2 

N l = 4n7 J ^ Ag » < 0, ’ X ') s 0) cos0* d 0' d \» 
0 ^rr/2 

if 

// 


Ag^ ^’» Xt ) c °s<2>' d (f* d X» 


Ag„ (<?\ x ') S«r) costf>’ d d X' 


(4) 
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The term can be calculated by numerical integration or by iterative inverse inter- 
polation (Reference 9). N is evaluated by a summation approximation over a 20°-by-20° 
square area centered at the evaluation point with data at a grid interval of l°-by-l° 
(mean free-air gravity anomalies are available for l°-by-l° square areas). Since the 
magnitude of the Stokes' function is known to be much smaller in area A than in A , 

N is taken to be negligible. 

O 

The accuracy of this calculation of N is estimated (Reference 9) to be ^2m on land and 
in coastal areas and * 7m in areas where surface gravity data is less dense. Figure 3 
shows a contour map for the North Atlantic area as a typical local geoid model. 

Generation of Power Spectra 

In general, three steps are involved in the generation of powe. spectra. These are the 
calculation of covariance functions, the application of a discrete cosine series transform 
to the covariance functions to generate "raw" power spectral density values, and the 
smoothing of these "raw" values to compensate for the effects of frequency analysis on 
a record of finite length. The following text treats these three steps in more detail. 

Blackman and Tukey (Reference 1G) give for autocovariance in one dimension at r 
arc degrees 



n-r 

C r " n^r ^ V N q+r ’ r = °» *• 2 * m ’ in ( me terc) 2 (5) 

q=l 


where m is the maximum number of lags or the largest correlation distance in terms 
of 1 arc degree intervals to be examined 


n is the total number of data points spaced at 1 arc degree intervals 


N 

q 


is the q geoid undulation, adjusted to zero mean on the set of n values 


For a two-dimensional array of data, such as in the detailed geoid models, the auto- 
covariance calculation is much more complicated, involving a sum of products of the 
undulation at a particular point p and the undulation at a point between r and 
r + 1 arc degrees distance from p. This sum is carried out over all points p and 
over all points within the (r, r + 1) distance criteria from each point p, with the 
proviso that no particular undulation product pair be represented in the summation more 
than once. This sum must be normalized by the number of undulation pair products 
included in the sum. The actual calculation is carried cut in a computer program and 
involves too many logical statements and tests to present here. However, the calcu- 
lation may be represented in a pseudo-mathematical form by 

c - — — — y y y N N t . (6) 

r p* c* P(r) Z-^ £—J Z-w q, s q's 

q=l s=l q’ s' 


where 


p is the number of rows in the data array 

c is the number of columns in the data array 


q', s' are the indices of any data point lying between r and r a 1 arc 
degrees distance from the point q, s 



indicates a summation over all points q' s' such that no duplication 

of products N N , , occurs in the total summation of equation (6). 

q, s q , s 
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Figure 3. DETAILED 


GRAVIMETRIC GEOID OF THE NORTH ATLANTIC 




P(r) 


is the number of data points q\s' included in the summation 
for a given value of r. 



This evaluation of the autocovariance is based on the ..ssumption that the autocovariance 
of geoid undulations is isotropic. 


Next a discrete finite cosine series transform is applied to the autocovariance alues 
to obtain a ’ raw 1 PSD function value, 


V = 2 
r 


m-1 

C + 2 T C cos + C cos rrr 
o q m m 

q=l 


m /cycle/arc degree (7) 


This raw estimate must be smoothed to compensate for discontinuities inherent in the 
digital discrete cosine transform process. The process known as "hamming" (Refer- 
ence 16) is used as follows, resulting in the smoothed power spectrum estimate, 

V = .54 V + .46 V, 

o o 1 

V ,23 Vl +,54 V* 23V M r - 1, 2, , m-1 (8) 

V = .46 V + .54 V 

m m-1 m 

Referring to equation (7) it is clear that a large maximum lag distance (large m), 
subject always to 

m < M 

where M - the smaller of p or c, 

will result in a large number of PSD estimates within the desired frequency band and 
therefore yield good frequency resolution. But PSD estimate accuracy decreases as 
the index r approaches m. The error E in the smoothed power spectrum (PS) is 
given (Reference 17) as a function of m and M as 

E - t/m/M* (9) 
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The use of large m values may also cause instability in the calculation of PSD v ilues 
as r approaches m. Instability in the PS plot manifests itself as periodic fluctuations 
in the PS as r approaches m. That is, the instability occurs at frequencies near the 
high frequency limit, which is the folding frequency 

V'iii < c y cle/de g-> < l0 > 

where h is the interval between data points, i. e., 1 arc degree. 

Exercising the tradeoff between resolution and accuracy or stability is largely a matter 
of cut-and-try in the choice of maximum lag value m, although Blackman and Tukey 
(Reference 16 ) recommend 

m^.lOM (11) 

After several PS plots for local geoids were generated, aliasing of the folding frequency 

was detected in the spectra. This effect is one of the pitfalls of sampling data at a 

regular interval and is manifest as artificial peaks in the power spectrum at rational 

fractions of the folding frequency, 

f f f 
c c c 

i.c., £ ’ 3 * 4 ’ etc. 

To remove this alias (to first order) power spectral density values were calculated, by 
the above method, using data whose power spectrum is known and which is similar to 
the geoid P.S. Aliased peaks in this calculated spectrum were easily identified and 
separated from the true spectrum. They were then normalized and subtracted from 
the calculated spectra for geoid heights. 

Kaula's Rule of Thumb 

This rule, derived from analysis of gravity anomalies (Reference 12 ), estimates the 
root-mean square value of the spherical harmonic geopotential coefficient of degree 
n from 
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( 12 ) 


a 

n 



si 


± 10 


-5 


n 


By expressing geoid undulation in spherical harmonics, the degree variance of geoid 
heights may be written (Reference 14) as 



fi 2 

- — (2n 1) a 

y n 


/ C . S l (meter ) 2 
( nin mu / 


(13) 


where n is the Earth’s gravitational constant and y is the average acceleration of 
gravity at the Earth’s surface. 


By substituting equation (12) into equation (13), it is seen that K aula’s rule-of-thumb 
implies 

_ 2 ux 10 10 / 2 n+l \ 0 

n { N / = ( 4 “ 1 (meter ) 2 (1 

The degree variance of N may be interpreted as the average amplitude squared of 
geoid undulation at a particular frequency, f n , 


where 


f 

n 


n 

360 arc degrees 


(15) 


it represents the value of the power spectra of geoid undulations at a particular fre- 
quency. If a continuum of frequencies is assumed than from equations 14 and 15 it is 
expected that the power spectra of geoid undulations V(f) would be approximately pro- 
portional to f -2 , 


V(f) 


K_ 

,3 


at high frequencies, if Kaula's rule-of-thumb remains valid. 


(16) 
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MODELING THE SPECTRA 


Figure 4 shows a pair of typical PSD plots, one for the area defined by the continental 
United States and one for the Northeast Pacific area. Both piois show good agreement 

_3 

with the f rule implied by Kaula's rule-of-thmcb. These plots also point up the 
general similarity between PSD plots of land and ocean geoid areas, a similarity that 
was observed for all nine geoid areas. This similarity is not expected to hold for 
frequencies much higher than 0. 5 cycles per arc degree due to the high amplitude 
undulations at short wavelengths expected in mountainous regions. The small peaks 
evident at the high frequency end of these plots are probably due to imperfect correction 
of the spectra for aliasing of the folding frequency. These plots were truncated at 0.4 
cycles per arc degree because of aliasing of higher frequency undulations at fre- 
quencies between 0.4 and 0. 5 cycles per arc degree. 

_3 

Figure 5 shows the consistency of the chosen f rule spectra model with all nine 

_3 

PSD plots. The placement of this f line for the derived spectra was adjusted to 
best approximate the PSD values of the two marine geoid areas represented, i.e. , 
the North Atlantic and the Northeast Pacific areas. The chosen model is represented 
by 


0.0525 2 


V = — 1 — - — m“/cycles/arc degree for f in cycles/arc degree ( 17 ) 
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m 2 /CYCLES/ARC DEGREE 


TYPICAL GEOID UNDULATION 
POWER SPECTRA 
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m 2 /CYCLES/ARC DEGREE 


DERIVED POWER SPECTRUM 
FOR GEOID PROFILES 



Figure 5 
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DISCUSSION OF RESULTS 


To assess the effect of the geoid spectra model on altimeter design and operation, 
assume an altimeter height resolution capability of ± lm. From Figure 5, a lm undu- 
lation (represented as lm^/cycle/ arc degree) is found to correspond to a frequency 
of 0.374 cycles/arc degree. Since 1 arc degree = 111 km., this frequency represents 
an undulation wavelength of 296 km. This represents, on average, the short wave- 
length limit for undulations that can be detected by an altimeter with 1 meter resolu- 
tion capability. However, actual resolution of this undulation also requires that the 
diameter of the surface area spot illuminated by the altimeter be much smaller than 296 
296 km., and the altimeter sampling rate must be such that at least two measurements 
be made in this distance. This implies a data rate of at least 2 measurements per 
minute at a typical satellite ground speed of 8 km./second. The GEOS-C altimeter, 
which has a planned 1 meter height resolution, should be able to detect and measure 
296 km. wavelength undulations in the North Atlantic since its effective illuminated 
spot diameter is about 12 km. and its sampling rate is at least 10 samples/second 
(Reference 18). 


CONCLUSIONS 

The power spectra of geoid undulations for the detailed geoid model are found to be 
not inconsistent with Kaula’s rule-of-thumb for the magnitude of harmonic coefficients 
cl different degree, following an inverse third power law with frequency. 

For an altimeter of lm precision operating over ocean areas, the geoid undulation 
power spectra imply requirements for at least a 2.0 measurement per minute data 
rate and a spot size of less than 50 km. diameter in order to resolve undulations of 
lm amplitude. 
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